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Quantum bits (qubits) are at the heart of quantum information processing schemes. Currently, 
solid-state qubits, and in particular the superconducting ones, seem to satisfy the requirements for 
being the building blocks of viable quantum computers, since they exhibit relatively long coherence 
times, extremely low dissipation, and scalability. The possibility of achieving quantum coherence 
in macroscopic circuits comprising Josephson junctions, envisioned by Legett in the 1980’s, was 
demonstrated for the first time in a charge qubit; since then, the exploitation of macroscopic quantum 
effects in low-capacitance Josephson junction circuits allowed for the realization of several kinds of 
superconducting qubits. Furthermore, coupling between qubits has been successfully achieved that 
was followed by the construction of multiple-qubit logic gates and the implementation of several 
algorithms. Here it is demonstrated that induced qubit lattice coherence as well as two remarkable 
quantum coherent optical phenomena, i.e., self-induced transparency and Dicke-type superradiance, 
may occur during light-pulse propagation in quantum metamaterials comprising superconducting 
charge qubits. The generated qubit lattice pulse forms a compound ’’quantum breather” that 
propagates in synchrony with the electromagnetic pulse. The experimental conhrmation of such 
effects in superconducting quantum metamaterials may open a new pathway to potentially powerful 
quantum computing. 
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Introduction 

Quantum simulation, that holds promises of solving particular problems exponentially faster than any classical 
computer, is a rapidly expanding field of research m- The information in quantum computers is stored in quantum 
bits or qubits, which have found several physical realizations; quantum simulators have been nowadays realized 
and/or proposed that employ trapped ions [1], ultracold quantum gases [S], photonic systems [0], quantum dots [7], 
and superconducting circuits DEli. Solid state devices, and in particular those relying on the Josephson effect 
[lOj . are gaining ground as preferable elementary units (qubits) of quantum simulators since they exhibit relatively 
long coherence times and extremely low dissipation m- Several variants of Josephson qubits that utilize either 
charge or flux or phase degrees of freedom have been proposed for implementing a working quantum computer; the 
recently anounced, commercially available quantum computer with more than 1000 ’ superconducting qubit CPU, 
known as D-Wave (the upgrade of D-Wave Two"^*'^ with 512 qubits CPU), is clearly a major advancement in 

this direction. A single superconducting charge qubit (SCQ) [T2] at milikelvin temperatures behaves effectively as 
an artificial two-level ’’atom” in which two states, the ground and the first excited ones, are coherently superposed 
by Josephson coupling. When coupled to an electromagnetic (EM) vector potential, a single SCQ does behave, 
with respect to the scattering of EM waves, as an atom in space. Indeed, a single-atom laser has been realized 
with an SCQ coupled to a transmission line resonator (cavity) [TJ]. Thus, it would be anticipated that a periodic 
structure of SCQs demonstrates the properties of a transparent material, at least in a particular frequency band. 
The idea of building materials comprising artificial ’’atoms” with engineered properties, i.e., metamaterials, and in 
particular superconducting ones El: is currently under active development. Superconducting quantum metamaterials 
(SCQMMs) comprising a large number of qubits could hopefully maintain quantum coherence for times long enough 
to reveal new, exotic collective properties. The first SCQMM that was only recently implemented comprises 20 flux 
qubits arranged in a double chain geometry m- Furthermore, lasing in the microwave range has been demonstrated 
theoretically to be triggered in an SCQMM initialized in an easily reachable factorized state m- 
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Results 


Superconducting Quantum Metamaterial Model 

Consider an infinite, one-dimensional (ID) periodic SCQ array placed in a transmission line (TL) consisting of two 
superconducting strips of infinite length [T71[TH] (Figure la and b); each SCQ, in the form of a tiny superconducting 
island, is connected to each bank of the TL by a Josephson junction (JJ). The control circuitry for each individual 
SCQ (Figure Ic), consisting of a gate voltage source Vg coupled to it through a gate capacitor Cg, allows for local 
control of the SCQMM by altering independently the state of each SCQ m- The SCQs exploit the nonlinearity 
of the Josephson effect and the large charging energy resulting from nanofabrication to create artificial mesoscopic 
two-level systems. A propagating EM field in the superconducting TL gives rise to nontrivial interactions between the 
SCQs, that are mediated by its photons [20) . Those interactions are of fundamental importance in quantum optics, 
quantum simulations, and quantum information processing, as well. In what follows, it is demonstrated theoretically 
that self-induced transparency m and Dicke-type superradiance (collective spontaneous emission) occur for weak 
EM fields in that SCQMM structure; the occurence of the former or the latter effect solely depends on the initial 
state of the SCQ subsystem. Most importantly, self-induced transparent (SIT) or superradiant (SRD) pulses induce 
quantum coherence effects in the qubit subsystem. In superradiance (resp. self-induced transparency), the initial 
conditions correspond to a state where the SCQs are all in their excited (resp. ground) state; an extended system 
exhibiting SRD or SIT effects is often called a coherent amlpifier or attenuator, respectively. These fundamental 
quantum coherent prosesses have been investigated extensively in connection to one- and two-photon resonant two- 
level systems. Superradiant effects have been actually observed recently in two-level systems formed by quantum dot 
arrays [ 23 ] and spin-orbit coupled Bose-Einstein condensates | 23 |; the latter system features the coupling between 
momentum states and the collective atomic spin which is analogous to that between the EM field and the atomic spin 
in the original Dicke model. These results suggest that quantum dots and the atoms in the Bose-Einstein condensate 
can radiatively interact over long distances. The experimental confirmation of SIT and SRD in extended SCQMM 
structures may open a new pathway to potentially powerful quantum computing. As a consequence of these effects, 
the value of the speed of either an SIT or SRD propagating pulse in a SCQMM structure can in principle be engineered 
through the SCQ parameters [25] . which is not possible in ordinary resonant media. From a technological viewpoint, 
an EM (light) pulse can be regarded as a ’’bit” of optical information; its slowing down, or even its complete halting 
for a certain time interval, may be used for data storage in a quantum computer. 

In the following, the essential building blocks of the SCQMM model are summarized in a self-contained manner, 
yet omitting unnecessary calculational details which are presented in the Supplementary Information. The energy per 
unit cell of the SCQMM structure lying along the x—direction, when coupled to an EM vector potential A = Az{x, t)z, 
can be readily written as [T7| |T8| 

H = ^{[^l-‘2cosipn] + [al + P^ian+i - anf] 

n 

-I-[2cos(p„(l - cosa„)]} , (1) 

in units of the Josephson energy Ej — $o-fc/(27rC'), with $ 0 ) Ic and C being the magnetic flux quantum, the critical 
current of the JJ, and the capacitance of the JJ, respectively. In equation Q, is the superconducting phase on the 
nth island, /3 = (87rdifj)“^/^($o/27r), with d being the separation between the electrodes of the superconducting TL, 
and the overdots denote differentiation with respect to the temporal variable t. Assuming EM fields with wavelengths 
X » £,d, with i being the distance between neighboring qubits, the EM potential is approximately constant within a 
unit cell, so that in the centre of the nth unit cell Az{x, t) ~ Az^„(t). In terms of the discretized EM potential A^^nit), 
the normalized gauge term is a„ = 2'KdAx^n/^o- The classical energy expression equation Q provides a minimal 
modelling approach for the system under consideration; the three angular brackets in that equation correspond to the 
energies of the SCQ subsystem, the EM field inside the TL electrodes, and their interaction, respectively. The latter 
results from the requirement for gauge-invariance of each Josephson phase. 


Second Quantization and Reduction to Maxwell-Bloch Equations 

The quantization of the SCQ subsystem requires the replacement of the classical variables and ipn by the 
corresponding quantum operators (pn and —i{d/d(pn), respectively. While the EM field is treated classically, the SCQs 
are regarded as two-level systems, so that only the two lowest energy states are retained; under these considerations. 


3 





■II-0 Vg 

Cg 

unit cell 


FIG. 1: Schematic drawing of a charge qubit superconducting quantum metamaterial (SCQMM). (a) The 

SCQMM comprising an infinite chain of identical charge qubits in a superconducting transmission line. Each qubit consists of a 
superconducting island that is connected to the electrodes of the transmission line through two Josephson junctions, formed in 
the regions of the dielectric layers (blue). The propagating electromagnetic vector potential pulse is also shown schematically 
out of scale, (b) The side view of the SCQMM in which the relevant geometrical parameters and the field orientations are 
indicated, (c) A unit cell of the superconducting quantum metamaterial which also shows the control circuitry of the charge 
qubit, consisting of a gate potential Vg applied to it through the gate capacitor Cg. 


the second-quantized Hamiltonian corresponding to equation Q is 

H = 'y ^ y ^ Ep[n)al^ pan,p + y ^ -|- /3 {oin+i — o-n) ] 

71 p n 

+4 (n)ajj sin , 

n p,p' 


( 2 ) 



















4 


where = 0 , 1 , i?o and Ei are the energy eigenvalues of the ground and the excited state, respectively, the operator 
al^ p {an,p) excites (de-excites) the nth SCQ from the ground to the excited (from the excited to the ground) state, 
and Vp^p' = / dtpE.p{tf) COB (p'E.p{tf) are the matrix elements of the effective SCQ-EM field interaction. The basis 
states Sp can be obtained by solving the single-SCQ Schrodinger equation — Ep + 2cos(p)Ep = 0. In 

general, each SCQ is in a superposition state of the form = J2p'^ri,p{t)al^ p\0). The substitution of into 
the Schrodinger equation with the second-quantized Hamiltonian equation (§ , and the introduction of the Bloch 
variables R^in) = 4'* RyM = Rzin) = |4'n.iP - l^n.oP, provides the 

re-formulation of the problem into the Maxwell-Bloch (MB) equations 

Rx{n) =—[A-I-81?sin^ i?p(n), 

Ry{n) = [A -I- 8i? sin^ ^ ] Rx{n) — 8p sin^ ^Rz{n), 

Rz{n) = -\-8p.B\v? ^Ryiri), 

that are nonlinearly coupled to the resulting EM vector potential equation 

an + + X [iJ-Rxin) + DRz{n)]} sina„ = /3^5a„, ( 6 ) 

where (5a„ = a„_i - 2a„ -h a„+i, D = (Vn - Eoo)/(2x): = (Voo + ^ii)/2, M = Eio/x = V^i/x, and A = d - Cq = 

(El — Eq)/x, with X = ^j/Ej- In the earlier equations, the overdots denote differentiation with respect to the 
normalized time t —> ujjt, in which uij = eIc/{T^C) is the Josephson frequency and e, h are the electron charge and 
the Planck’s constant devided by 27r, respectively. 


( 3 ) 

( 4 ) 

( 5 ) 


Approximations and Analytical Solutions 


For weak EM fields, the approximation sina„ ~ a„ can be safely used. Then, by taking the continuum limit 
cin{t) — >■ a{x, t) and Ri{n-, t) —>■ Ri{x] t) {i = x, y, z) of equations (3][5 ) and Q, a set of simplified, yet still nonlinearly 


coupled equations is obtained, similar to those encountered in two-photon SIT in resonant media m- Further 
simplification can be achieved with the slowly varying envelope approximation (SVEA) by making for the EM vector 
potential the ansatz a(x,t) = £{x,t) cos'^(x,t), where '^{x,t) = kx — ujt 4>{x,t) and £{x,t), (l>{x,t) are the slowly 
varying pulse envelope and phase, respectively, with w and k = being the frequency of the carrier wave 

of the EM pulse and its wavenumber in the superconducting TL, respectively. In the absence of the SCQ chain the 
EM pulse is ’’free” to propagate in the TL with speed /3. At the same time, equations (|3][^ for the Bloch vector 
components are transformed according to Rx = Tx cos ( 24 ') -|- Xy sin(24'), Ry = Xy cos(24') — sin(24'), and Rz = Xz- 

Then, collecting the coefficients of sin dt and cos dt while neglecting the rapidly varying terms, and averaging over 
the phase 4', results in a set of truncated equations (see Supplementary Information). Further manipulation of the 
resulting equations and the enforcement of the two-photon xesonance condition A = 2w, results in 


£-\-C£x = 
(j) -\- C(j)x = 


A* 

-X^sry, 

2D 

-X^Tz, 


( 7 ) 

( 8 ) 


where c = j3"^k/uj = 2j3'^k/IX, and the truncated MB equations 


Xx = -2D£^Xy, Xy = -\-2 D£^Xx - Xz = +^ry, ( 9 ) 

which obey the conservation law -|- = 1 . In equations the n—dependence of the Xi {i = x,y,z) is 

suppressed, in accordance with common practices in quantum optics. 

The Xi can be written in terms of new Bloch vector components Si using the unitary transformation Xx = Sx cos 4> — 
E^sin^, Xy = Sy, and Xz = S' 2 C 0 s 4 > -|- S'a;Sin$, where $ is a constant angle to be determined. Using a procedure 
similar to that for obtaining the Xi, we get Sx = 0, Sy = — ^We'^Sz, and Sz = -\-^W£‘^Sy, where LE = a/ (AD)^ -f 
and tan4> = 7 = AD/fi. The combined system of the equations for the Si and equations Q-(|^ admits exact solutions 
of the form e = e(r = t — x/v) and Si = Siir = t — x/v), where v is the pulse speed. For the slowly varying pulse 
envelop, we obtain 

I. 

2 



e{T) = So 


( 10 ) 
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FIG. 2: The velocity-amplitude relation in two-photon superradiant (TPSRD, amplifying) and two-photon 
self-induced transparent (TPSIT, absorbing) superconducting quantum metamaterials (SCQMMs) & quantum 
coherent pulse profiles. In all subfigures, the pulse velocity v in units of /3 as a function of the electromagnetic vector potential 
pulse amplitude eq is plotted and compared with the corresponding curves for ordinary (atomic) amplifying and absorbing 
media (brown- and green-dotted curves, respectively). The horizontal magenta-solid (resp. black-solid) lines indicate the 
limiting velocity in ordinary amplifying and absorbing media, v/P = 1 (resp. amplifying and absorbing SCQMMs, v = c < /3). 
(a) Vbo = Vii = 1, Voi = V\o = 0.8, x = 1/5, l?i — = 3 (7 = 0 and fl/w = 0.3). Left Inset: The electromagnetic vector 

potential pulse envelop (e/ea/)^ and the population inversion function Rz{n) profiles as a function of the slow variable {t/tm 
in a frame of reference that is moving with velocity v, for TPSIT (absorbing) SCQMMs. Right Inset: Same as in the left inset 
for TPSRD (amplifying) SCQMMs. (b) Vbo = 0.6, Vii = 1.4, Vbi = Vbo = 0.8, x = 1/5, i?i — So = 3 (7 = 2 and Q,/ui = 0.3). 
Left Inset: The electromagnetic vector potential pulse envelop (ejeM)^ and the population inversion function Rz{n) profiles as 
a function of the slow variable (j jrM in a frame of reference that is moving with velocity u, for TPSIT (absorbing) SCQMMs in 
the presence of relatively strong decoherence (7 = 2). Right Inset: Same as in the left inset for TPSRD (amplifying) SCQMMs. 
(c) Vbo = Vii = 3, Vbi = Vio = 0 . 8 , x = 1/5, Si — So = 3 (7 = 0 and D/tu = 0.52). (d) Vbo = 3, Vii = 3.8, Vbi = Vio = 0 . 8 , 
X = 1/5, Si — So = 3 (7 = 2 and fl/o; = 0.52). The effect of non-zero decoherence (7 yf 0) become apparent by direct 
comparison of a with b and c with d. The pulse velocity v in SCQMMs saturates with increasing eo to Vm/ P, that can be 
significantly lower than that achieved in ordinary TPSIT and TPSRD media, i.e., p. The parameters of the SCQMM can be 
engineered to slow down the pulse velocity v at the desired level for high enough amplitudes Eo. Note that v is also the velocity 
of the coherent qubit pulse. 


where Eq = ■>/ /lo)[v/{ c — u)] is the pulse amplitude and Tp = {x{u^./ijj)[v/{c — v)]} ^ its duration, with a = 
m/vf = 1/7/r + 7 ^. The decoherence factor 7 can be expressed as a function of the matrix elements of the SCQ-EM 
field interaction, V/j, as 7 = 2(Vii — Voo)/hio that can be calculated when the latter are known. Such Lorentzian 
propagating pulses have been obtained before in two-photon resonant media |281I29) : however, SIT in quantum systems 
has only been demonstrated in one-photon (absorbing) frequency gap media, in which solitonic pulses can propagate 
without dissipation [30] . The corresponding solution for the population inversion, reads 


R.{t) = ± 



( 11 ) 


where Em = 2-^(l/a;)[u/(c — u)], and the plus (minus) sign corresponds to absorbing (amplifying) SCQMMs; these 
are specified through the initial conditions as Rz{—oo) = —1, e(— 00 ) = 0 and i ?^(—00 = -1-1), e(— 00 ) = 0 for 
absorbing and amplifying SCQMMs, respectively (with Rx{—oo) = Ry{—oo) = 0 in both cases). The requirement for 
the wavenumber k being real, leads to the SCQ parameter-dependent condition 2x^(Cii -I- Vbo) < {Ei — EqY loi' pulse 
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propagation in the SCQMM. Thus, beyond the obtained two-photon SIT or SRD, the propagating EM pulse plays a 
key role in the interaction processes in the qubit subsystem: it leads to collective behavior of the ensemble of SCQs 
in the form of quantum coherent probability pulses; such pulses are illustrated here through the population inversion 

R.. 

The corresponding velocity-amplitude relation of the propagating pulse reads 


1 


Loe, 


0 J 


( 12 ) 


Equation (12) can be also written as a velocity-duration expression, since the pulse amplitude and its duration are 


related through s^Tp = 4/(xlE)- The duration of SRD pulses cannot exceed the limiting value of tm = w(c — 
v)/{xfiv). From equation ( [l^ , the existence of a critical velocity c, defined earlier, can be immediately identified; 
that velocity sets an upper ^wer) bound on the pulse velocity in absorbing (amplifying) SCQMM structures. Thus, 
in absorbing (amplifying) SCQMM structures, pulses of higher intensity propagate faster (slower). That limiting 
velocity is generally lower than the corresponding one for two-photon SIT or SRD in ordinary media, /3, which here 
coincides with the speed of the ’’free” pulse in the TL (Figure 2). As can be inferred from Figure 2, the increase 
of decoherence through 7 makes the velocity to saturate at its limiting value c at lower amplitudes e; that velocity 
can be reduced further with increasing the ratio of the TL to the pulse carrier wave frequency Q/uj through proper 
parameter engineering. Moreover, effective control of v in SCQMMs could in principle be achieved by an external field 
m or by real time tuning of the qubit parameters. That ability to control the flow of ’’optical”, in the broad sense, 
information may have technological relevance to quantum computing [25]. Note that total inversion, i.e. excitation or 
de-excitation of all qubits during pulse propagation is possible only if 7 = 0 , i.e., for V()o = Vu; otherwise (Vqo < Vn) 
the energy levels of the qubit states are Stark-shifted, violating thus the resonance condition. Typical analytical 
profiles for the EM vector potential pulse e(r) and the population inversions Rz(j) both for absorbing and amplifying 
SCQMMs are shown in the insets of Figure 2. The maximum of e(r) reduces considerably with increasing 7 , while at 
the same time the maximum (minimum) of Rz decreases (increases) at the same rate. 

The system of equations 0 -® and ([^ can be reduced to a single equation using the parametrization = 
i?o 7 CT^[l — cosd], ry = —Rocsind, and = i?o{l — cr^[l — cos0]}, of the Bloch vector components. Then, a relation 
between the Bloch angle 9 = 9{x,t) and the slow amplitude e can be easily obtained, that leads straightforwardly 
to the equation 0 + c6x = —i?oX^^cos0. Time integration of that equation yields || = RoX'^i^ ~ cos 9), that 
conforms with the famous area theorem: pulses with special values of ’’area” 9{x) = 27rn conserve that value during 
propagation. 

Here we concentrate on the interaction of the SCQs with the EM wave and we are not concerned with decoherence 
effects in the SCQs due to dephasing and energy relaxation. This is clearly an idealization which is justified as long as 
the coherence time exceeds the wave propagation time across a relatively large number of unit cell periods. In a recent 
experiment [ 2 S], a charge qubit coupled to a strip line had a dephasing time in excess of 200 ns, i.e., a dephasing 
rate of 5 MHz, and a photon loss rate from the cavity of 0.57 MHz. Those frequencies are very small compared 
with the transition frequency of the considered SCQs which is of the order of the Josephson energy (i.e., a few GHz) 
[nittHj- Therefore, we have neglected such decoherence effects in the present work. The decoherence factor 7 , which 
in Figures 2b and 2d has been chosen according to the parameter values in m, is not related to either dephasing or 
energy relaxation. That factor attains a non-zero value whenever the matrix elements of the effective SCQ-EM field 
interaction, Vu and Vboj are not equal. 


Numerical Simulations 

In order to confirm numerically the obtained results, the equations 0-0 and 0 are integrated in time using a 
fourth order Runge-Kutta algorithm with constant time-step. For pulse propagation in absorbing SCQMMs, all the 
qubits are initially set to their ground state while the vector potential pulse assumes its analytical form for the given 
set of parameters. A very fine time-step and very large qubit arrays are used to diminish the energy and/or probability 
loss and the effects of the boundaries during propagation, respectively. The subsequent temporal evolution in two- 
photon SIT SCQMM, as can be seen in Figures 3a and 3b, in which several snapshots of the population inversion 
Rz{n;t) and the vector potential pulses an{t), respectively, are shown, reveals that the latter are indeed capable of 
inducing quantum coherent effects in the qubit subsystem in the form of population inversion pulses! In Figure 3a, 
the amplitude of the Rz{n;t) pulse gradually grow to the expected maximum around unity in approximately 60 time 
units, and they continue its course almost coherently (although with fluctuating amplitude) for about 160 more time 
units, during which they move at the same speed as the vector potential pulse (Figure 3b). However, due to the 
inherent discreteness in the qubit subsystem and the lack of inter-qubit coupling, the i?^(n;t) pulse splits at certain 
instants leaving behind small ’’probability bumps” that get pinned at particular qubits. After the end of the almost 
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FIG. 3: Numerical validation of the analytical expressions for two-photon self-induced transparent (TPSIT) 
and superradiant (TPSRD) propagating pulses, (a) Snapshots of the population inversion pulse Rz{n-,t), excited 
by the induced quantum coherence in the qubit subsystem by the electromagnetic vector potential pulse, in the absence of 
decoherence (7 = 0); the pulse propagates to the right (time increases downwards) in TPSIT (absorbing) superconducting 
quantum metamaterials (SCQMMs). The snapshots are taken at intervals of 20 time-units starting at t = 20 and they are 
displaced vertically to avoid overlapping (blue pulses). The corresponding pulses from the analytical expression equation © at 
the same time-instants are shown in red. (b) Snapshots for the corresponding evolution of the electromagnetic vector potential 
pulse a„{t), that exhibits significant broadening as time passes by; the numerical and analytical pulses are shown in blue and 
red color, respectively, (c) The same as in a in TPSRD (amplifying) superconducting quantum metamaterials. The resulting 
propagation is not as simple as expected from the theoretical analysis; instead of a population inversion pulse, it is observed 
a rather kink-like front propagating to the the right (blue) with a velocity considerably less than that predicted analytically 
for the pulse, which analytical form is shown in red. (d) The same as in b in TPSRD (amplifying) superconducting quantum 
metamaterials. The velocity of the an{t) pulse (blue) is the same as that of the propagating population inversion front, Rz{n-, t); 
however, it exhibits less broadening with time in comparison with the corresponding numerical a„(t) pulse in b. The predicted 
analytical form is shown in red. Parameter values: x = 1 / 5 , /3 = 6 , Vbo = Fii = 1, Vbi = Vio = 0.8, Ei— Eq = 3, and vlc = 0.7 
(for a and b); v/c = 1.25 (for c and d). 


coherent propagation regime, the Rzin; t) pulse broadens and slows-down until it stops completely. At the same time, 
the width of the a„(t) pulse increases in the course of time due to discreteness-induced dispersion. A comparison 
with the corresponding analytical expressions reveals fair agreement during the almost coherent propagation regime, 
although both the Rz{n\t) and a„(f) pulses travel slightly faster than expected from the analytical predictions. The 
temporal variable here is normalized to the inverse of the Josephson frequency wj which for typical parameter values 
is of the order of a few GHz m- Then, the almost coherent induced pulse regime in the particular case shown in 
Figure 3 lasts for ~ 160 x 10“® s, or ~ 160 ns, which is of the same order as the reported decoherence time for a 
charge qubit in [5S] (i.e., 200 ns). 

The situation seems to be different, however, in the case of two-photon SRD pulses, as can be observed in the 
snapshots shown in Figures 3c and 3d for Rz{n] t) and an{t), respectively. Here, the lack of the inter-qubit interaction 
is crucial, since the SCQs that make a transition from the excited to their ground state as the peak of the a„(t) 
pulse passes by their location, cannot return to their excited states after the a„(t) pulse has gone away. It seems, 
thus, that the a„(t) pulse creates a type of a kink-like front that propagates at the same velocity. It should be noted 
that the common velocity of the Rz(n-,t) and a„(t) pulses is considerably different (i.e., smaller) than the analytically 
predicted one, as it can be inferred by inspection of Figures 3c and 3d. Even more complicated behavioral patterns of 
two-photon SRD propagating pulses and the effect of non-zero decoherence factor are discussed in the Supplementary 
Information. 














































































































Conclusion 


An SCQMM comprising SCQs loaded periodically on a superconducting TL has been investigated theoretically 
using a minimalistic one-dimensional model following a semiclassical approach. While the SCQs are regarded as 
two-level quantum systems, the EM field is treated classically. Through analytical techniques it is demonstrated that 
the system allows self-induced transparent and superradiant pulse propagation given that a particular constraint is 
fulfilled. Most importanty, it is demonstrated that the propagating EM pulses may induce quantum coherent pop¬ 
ulation inversion pulses in the SCQMM. Numerical simulation of the semiclassical equations confirms the excitation 
of population inversion pulses with significant coherence time in absorbing media. The situation is slightly differ¬ 
ent in amplifying media, in which the numerically obtained, induced population inversion excitations are kink-like 
propagating structures (although more complex behaviors discussed in the Supplementary Information also appear). 
Moreover, the limiting pulse velocity in both amplifying and absorbing SCQMMs is lower than the corresponding one 
in two-photon resonant amplifying and absorbing ordinanary (atomic) media. That limiting velocity in SCQMMs can 
in principle be engineered through the SCQ parameters. 
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I. HAMILTONIAN FUNCTION AND QUANTIZATION OF THE QUBIT SUBSYSTEM 

Consider an infinite number of superconducting charge qubits in a transmission line (TL) that consists of two 
superconducting plates separated by distance d, while the center-to-center distance between qubits, £, is of the same 
order of magnitude. The charge qubit is of the form of a mesoscopic superconducting island which is connected to each 
electrode of the TL with a Josephson junction (JJ). Assume that an electromagnetic (EM) wave corresponding to a 
vector potential A = A^ix, t)z propagates along the superconducting TL, in a direction parallel to the superconducting 
electrodes and perpendicular to the direction of the EM wave propagation. A minimalistic description of that system 
constists of writing the Hamiltonian of the compound qubit array - EM field system with the geometry shown in Fig. 
1 of the paper, as 


= X! - 2cosa„ cos^pn + al+ /3^(a„+i - Unf) , 

n 


(13) 


where ipn is the superconducting phase on n—th island, /3^ = {8TTdEj)~^ ($o/(27r))^, with d being the separation 
between the superconducting electrodes of the TL, and the overdots denote derivation with respect to the temporal 
variable t. In what follows, a basic assumption is that the wavelength A of the EM field is much larger than the other 
length scales determined by the size of the unit cell of the qubit metamaterial (A >> i,d), so that the vector potential 
component A^{x,t) can be regarded to be approximatelly constant within a unit cell, Az{x,t) — A^^nit). Then, 


the discretized and normalized EM vector potential at the n—th unit cell which appears in Eq. (13) reads a„(t) = 
(27rd/<I>o)^n(f)- The Hamiltonian function Eq. (13) is given in units of the Josephson energy Ej = (<i)o/c)/(27rC'), 
where Ic and C is the critical current and capacitance, respectively, of the JJs, and = h/(2e) is the flux quantum, 
with h and e being the Planck’s constant and the electron charge, respectively. That Hamiltonian can be written in 
a more transparent form by adding and subtracting 2 cos (()„ and subsequently rearranging to get 


H = H, 


QUB 


H 


EMF ' 




(14) 


where the qubit subsystem energy Hqjjb, the EM field energy Hemf, and their interaction energy Hint, take 
respectively the form 


HquB = - 2 cos (fin}, 


Hemf = 


+/3‘^{an+l - anf}, H,nt = y^{2cOStpn(l - COS On)}. (15) 


In the following, the EM field is treated classically, while the qubits are regarded as two-level systems. The latter 
approximation is particularly well suited in the case of resonant qubit - EM field interaction adopted here. 

The quantization of the qubit subsystem can be formally performed by replacing the classical variables pn and 
Pn by the quantum operators pn and pn —t respectively, in the Hamiltonian Hqub since we are dealing 

with a large number of Cooper pairs. The exact energy spectrum Ep(n) and the corresponding wavefunctions 'E.p{n) 
of the nth qubit may then be obtained by mapping the Schrodinger equation with the single-qubit Hamiltonian 
Hsq = Pn ~ cos Pn , onto the Mathieu equation 


52 

F 2 T £dp n 2 cos Pn 

dpi 


-‘P,n 


= 0. 


(16) 


The second quantization of the qubit subsystem proceeds by rewritting the single-qubit Hamiltonian as 


H. 


sq 


Hsq = -f 2cos(^^4'(v?), 


(17) 


where and 'F are field operators. Note that the subscript ”n” in Eq. (17) has been dropped since the qubits are 
identical. Using the expansion ^(yi) = J^p^p'^pip)^ where the operators (op) create (annihilate) qubit excitations 
of energy Ep, the Hamiltonian Eq. (Il7|) is transformed into 


H,g = 


^ ^ EpOpOp. 
p=0,l,... 


(18) 
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We hereafter restrict Hgq to the Hilbert subspace of its two lowest levels, i.e., those with p = 0,1, so that in second 
quantized form the Hamiltonian Eq. ( |14[ ) reads 

H = EE E{ '^n +/3^(Q^n+l ~ Q^n)^}) (19) 

n p PiP' 


where p,p' = 0,1 and the lp,p' (n) = Vpi^p{n) that represent the matrix elements of the nth qubit - EM field interaction, 
are given by 


^PiP'i'^) — J' ^‘•Pn^pi^n') ‘•Pn^p,n{^n) ■ 


( 20 ) 


In the reduced state space, in which a single qubit can be either in the ground {p = 0) or in the excited {p = 1) state, 
the normalization condition al^ pUn^p = 1 holds for each qubit in the metamaterial. The reduced Hamiltonian Eq. 
(19) is also written in units of Ej. 


II. MAXWELL-BLOCH FORMULATION OF THE DYNAMIC EQUATIONS 


In accordance with the adopted semiclassical approach, the time-dependent Schrodinger equation 

= ( 21 ) 

in which H is the Hamiltonian from Eq. (jl^ in physical units, i.e., H = HEj, is employed for the description of the 
qubit subsystem. Generally, the state of each qubit is a superposition of the form 

( 22 ) 

p 

for which it can be easily shown that the coefficients satisfy the following normalization conditions 


E 


^n,p{t)? = 1 , 


Y,\^nAt)\^=N, 


(23) 


where in the second Eqs. (23) a finite A—qubit subsystem has been implied. Assuming that the pulse power is not 
very strong, substantial simplification of the dynamic equations for the qubit subsystem can be achieved using the 
approximation [1 — cos(a„)] ~ (l/2)a^ in the interaction Hamiltonian Substitution of |d>) = from Eq. 


(22) into the Schrodinger equation (21), along with derivation of the classical Hamilton’s equation for the normalized 
EM vector potential yields 


— ^p^n,p ~b ^ (^)5 

dn ^ “h CXji—l -f ^ ^ — 0, 

p,p' 


(24) 

(25) 


where x = ‘huij/Ej. In Eqs. (24) and (25), the temporal variable is renormalized according to t 
of that, the dimensionless energy of the qubit excitations has to be redefined according to Ep —>■ e 


—>■ ojjt. Because 
P = Ep/x, which 


explains the presence of the dimensionless factor l/y in Eq. (24). 

The evolution Eqs. (24) and (25) can be rewritten in terms of the n—dependent Bloch vector components through 
the transformation 


Rz{n) = - l^'n.ol 


Ry[n) = i(4':,o^n.i - «':.i«'n.o), Rx{n) = 


(26) 


in which the variables i?j (i = x,y,z) apply to each single qubit. Thus, from the point of view of a macroscopic 
system, represent the population densities Qi = By transforming Eqs. (24) and (25) according to 

Eq. (26) we get 

Rx{n) = -(A + 2Da^^)Ry{n), Ry{n) = +(A + 2Da^)Rx{n) - Rz{n) = (27) 

“1“ X[^ ~ ^ (o^n—1 ^n+l); 














where 
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_ (Vii — Vbo) ^2 _ (^11 + K)o) 

“ 2 x ’ " 2 ^ 


^10 A (El — Eq) 

^JL= -, A = ei-eo= -. 

X X 


(29) 


By taking the continuous limit of Eqs. (27) and (28) using the relations a„ —>■ a{x,t) and a„±i a ± ax + ^axx, we 
obtain 


Rx = -(A + 2Da^)Ry, Ry = +(A + 2Da^)Rx - 2^a^R^ 


R^ — - 1 - 2 ^ 0 ; Ryi 


a - 13‘^axx + = -xiDR^ + iJ.Rx)a, 


(30) 

(31) 


where Rx, Ry, Rz, and a are functions of the spatial variable x and normalized temporal variable t, while the overdots 
denote partial derivation with respect to the latter. It can be easily verified that the Bloch equations Eqs. (30) possess 
the dynamic invariant i?| = 1, whose constant value is derived through the normalization conditions. 


III. SLOWLY VARYING ENVELOPE APPROXIMATION (SVEA) 

The Slowly Varying Envelope Approximation (SVEA) relies on the assumption that the envelop of a travelling 
pulse in a nonlinear medium varies slowly in both time and space compared with the period of the carrier wave, that 
makes possible the introduction of slow and fast variables. Assuming that in the following, we can approximate the 
EM vector potential function by 


a{x,t) = e{x,t) costp{x,t), 


(32) 


where ip{x, t) = kx — ujt + (j)(x, t), with k and w being the wavenumber and frequency of the carrier wave, respectively, 
that are connected through the dispersion relation that is obtained at a later stage. The functions £{x,t) and 4>{x,t) 
in Eq. (32) are the slowly varying envelop and phase, respectively. Using fast and slow variables, the x and y Bloch 
vector components, Rx{n) and Ry{n), can be expressed as a function of new, in-phase and out-of-phase Bloch vector 
components Vx and Xy as 


Rx = Xx cos(27/i) -I- Ty sm{2'ijj), Ry = Xy cos(2^/i) — Xx sin(2i/'). 


Rz = Tz 


(33) 


A. Derivation of the dynamic equations for the slowly varying envelop and phase of the electromagnetic 

vector potential 


From Eqs. (32) and (33), the second temporal and spatial derivative of a(x,t) can be approximated by 
a « 2a;e sin ijj + {2ijj(f> — uj'^)£ cos ip, axx ~ —2-kSx sin ip — {2k(px — k^)£ cos ip, 


(34) 


in which the rapidly varying terms of the form e, £xx, (p"^, 4‘xx, <P, 4>x^x, etc., have been neglected. Substitution of 
Eqs. (|M|) and (33), into Eq. (31) gives 


2(u)£ + kp‘^£x) sinip + [2{(puj + kcpx) — ui'^ + + /3^fc^]ecos'(/' = —x{^U + k'ixx cos{2ip) + Xy sin(2^)]}ecos'0- (35) 


Equating the coefficients of sin ip and cos ip in the earlier equation yields 

we -I- kjS’^Ex = —xk-XyS cos^ ip, 

and 

2(0w -I- k(px) — {w^ — — P^k^} = —x[Dxz + yxx cos(2V')]. 


(36) 


(37) 


In order to obtain the dispersion relation w = uj(k), we require the vanishing of the expression in the curly brackets 
in Eq. (37), which yields for the wavevector of the EM radiation (i.e., the pulse) in the TL the expression 


k = ± 


Vw 2 - 

■ 


(38) 


Thus, the EM wave can propagate through the superconductin g q uantu m m etamaterial only when its frequency 
exceeds a critical one, Wc = U = V^O^ocT+^^TiV^- Finally, Eqs. (36) and (37) are averaged in time over the period 
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T = 27r/a; (i.e., the fast time-scale) of the phase Due to the simple time-dependence of tp{x,t) that has been 

assumed earlier in the framework of SVEA, that averaging actually requires the calculation of integrals of the form 


1 

{J^{smf{tf}),cosg{tpj)) = — J J^{smf{tp),cosg{ip))dtp. 


This procedure, when it is applied to the two evolution equations Eqs. 
for slow amplitude and phase 


(36) and (37) gives the truncated equations 


JL 

~2uj' 


> + C(j)x= -X — Rz, 

UJ 


where c = P'^k/oj is a, critical velocity. 


(39) 


B. Derivation of dynamic equations for the transformed Bloch vector components 


Substituting Eqs. ( |32[ ) and ( |33[ ) into Eqs. ( |30[ ) for the original Bloch vector components, we get 

(f-a; -I- 2'tpry) cos(2i/)) -I- {vy — 2iproc) sin(2^) = — (A -|- 2De^ cos^ '^)[i"y cos(2'!/') — Tx sin(2'!/))] 


(40) 


{r-y — 2'tprx) cos(2'0) — (r^, -I- 2'ipry) sin(2?/;) = (A -|- 2Ds^ cos^ '4’)[rx cos(2'(/i) -|- Vy sin(2'!/;)] — 2/ie^ cos^ ipr^ (41) 


Vz = 2g,e^ cos^ '(/'[’’y cos(2'0) — sin(2'0)]. 


(42) 


In order to derive the dynamic equations for the new Bloch vector components r^, we first multiply Eqs. (40) and 
(41) by cos(2'0) and sin(2-0), respectively, and then subtract one equation from the other. Thus we get 


Tx -\- ‘2'ipry = — (A -I- 2De^ cos^ cos^ if) cos(2'0)r^. 


Similarly, by multiplying Eqs. (40) and ([4l|) by sin(2'(/i) and cos{2ip), respectively, and by adding them, we get 


Vy — 2'tprx = (A -I- 2De^ cos^ 'p)rx — ‘2gs^ cos^ 'p sm{2p)rz- 


(43) 


(44) 


Finally, an average over the phase p is performed onto Eqs. (42 )-(|44| using the (easy to prove) relations 
(cos^')/'cos(2'(/))) = 1/4 and (cos^ sin(2^)) = 0. That procedure, and using also the relation ip = p — oj, yields 
the truncated Bloch equations 




— —(<5 -l- 2p -|- Ds )ry — —£ 


fy — +(<5 + 2,(f> -\- DP^^Tx 


d 2 

rz = 2 ^ ry, 


(45) 


where 5 = A — 2a;. Note that Eqs. (45) possess a dynamic invariant that has a form similar to that of the original 
Bloch equations, Eqs. (30), i.e.. 


2 I 2 I 2 1 

rx+ry + T^ = 1 . 


(46) 


The truncated Bloch equations Eqs. (45), along with Eqs. (39) for the slowly varying amplitude and phase of the 


EM vector potential pulse, constitute a closed system with five unknown functions which describe the approximate 
dynamics of the superconducting quantum metamaterial that are determined in the next section. 


IV. EXACT INTEGRATION OF THE TRUNCATED EQUATIONS 


The combination of Eqs. (39) and the third equation of Eqs. (45) provides a relation between the slow amplitude 


and the phase of the EM vector potential pulse. The first of Eqs. (39) is multiplied by e and written as 


^ A 


£^{x,t) = -x-£^{x,t)r, 

UJ 




(47) 
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(451. Thus we get 


while the derivative of the second of Eqs. (39) is taken with respect to time and then = rz is replaced from Eqs. 


^ d_ 

dt^'^dx 


= -x^s‘^{x,t)ry. 


(48) 


From Eqs. (47) and (48), and by taking into account the independence of the slow temporal and spatial variables, we 
get 


2(j>(x, t) = De^{x, t) + const., 


(49) 


where the constant of integration can be set equal to zero. Using Eq. (49), the truncated Bloch equations Eqs. (45) 
can be written as 


= -{5 + 2De^)ry, fy = +(<5 + 2D£^)r,z - ^s^rz, 


M 2 

= 2 ^ ^y 


' y, I y , ,-j, 2 

The latter can be written in a simpler form using the unitary transformation 

= Sx cos $ — S';; sin $, Ty = Sy, = Sz COS 4) + Sa; sin 4) 


(50) 


(51) 


where 4> is the constant transformation angle that is going to be determined. The truncated Bloch equations for the 
Ti, i = x,y, z, can be written in terms of the new Bloch vector components Si, using a procedure similar to that used 
in Subsection III.B to obtain Eqs. (45). Sustituting Eqs. into Eqs. ([^, we get 


Sx cos 4) — Sz sini^ = —{5 + 2De^)Sy, 


(52) 


Sy = 


[5 + 2De^) cos 4> — — sin 4? 


Sx- 


(<5 + 2D£^) sin 4) + —cos 4) 


5. 


(53) 


Sx sin 4) + Sz cos 4> = —^e^Sy. 


(54) 


Multiplying Eqs. (52) and (54) by cos4> and sin 4), respectively, and then adding them together, we get 

Sx = <e^ 


1 


- y. sin 4' — 2D cos 4> 


— S cos 4) )■ Sy. 


(55) 


Similarly, multiplying Eqs. (52) and (54) by sin4> and cos4>, respectively, and then subtracting one equation from 
the other, we get 


Sz = <e 


-fi cos dJ + 2D sin 4> 


— 5sin4) Sy. 


(56) 


The transformation angle $ can now be selected so that the resulting equations are simplified as much as possible. 
We thus define 


tan 4> = 7 = 


4D 

7 


(57) 


so that cos4> = ±cr and sin4> = icy where a = 1 /a/1 + y^. The choice of the sign is irrelevant and here we 
pick positive sign for both functions. Us in g th e abo ve v alue for the transformation angle and the definitions W = 
1 /( 474)2 _|_ ^2 Eqs. (55), (53), and (56) obtain their final form 


Sx - “t" ySy 


Sy - ySx 


1 9 

777 - Ws 

' ’ 2 


Sz, Sz = 


-r]l+\w£^ 


Sy. 


(58) 


In order to investigate the possibility for ’’coherent propagation” of an electromagnetic potential pulse, we consider 
resonance conditions. In that case, rj — 0, and Eqs. (58) become 


= 0 , 


Sy = -\we^Sz, Sz = +\w£^Sy. 


( 59 ) 
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Combining the second and third Eq. (591 and integrating, we obtain the ’’resonant” conservation law Sy+Sl = const.. 
Assuming that all the qubits are in the ground state at t = —oo, we have the initial conditions rx{t = —oo) = 
Tyit = —oo) = 0 and Czit = —oo) = —1 which are transformed into Sx(t = —oo) = Sy{t = —oo) = 0, and 

Szit = —oo) = —a through the transformation Eq. (51). Applying the initial conditions to the resonant conservation 
law, we get 


Si = nC 


(60) 


In the following, we seek solution of the form e = £(t = t — x/v) and Si = Si(T = t — x/v), with i = x,y,z. By 
changing the variables in the first of Eqs. (39), with Vy being replaced by 5^, we get after rearramgement 


Er fJ- V ^ 

T=x-—S.. 


~2uj c — V 


Then, combining Eq. (41) with the third Eq. (59) and integrating, we get 


(61) 


(62) 


where the conditions £(—oo) and Szi—oo) = —a were used. The system of Eqs. (60)-(62) for e, Sy, and Sz can be 
integrated exactly; the variables Sy and Sz can be eliminated in favour of £ to give £,- = Xe^ Va + be '^, in which the 
constants are defined as a = 2a/ k, b = A = simplify the notation. The equation for 

£ can be readily integrated 


de 


£2\/a + be^ 


= X 


, Va + be'^ . 

dr ^ - = A(t - To), 


ae 


(63) 


where we have set Eq = e(t = tq) = \/2aK to eliminate the boundary term resulting from the integral over e. Solving 
Eq. (63) for £, we finally get a Lorentzian-like slowly varying amplitude 


e(t) = 


£o 


]/l+Tp^(T-To)^ 

where £o = ^/-a/b = ^f2aK and ■ 


(64) 


V. DETAILS OF THE NUMERICAL SIMULATIONS AND THE ROLE OF DECOHERENCE 

The numerical simulations in order to confirm (or not) the existence of the quantum coherent effects and the 
quantum coherence induced in the medium, i.e., the qubit transmission line, which are predicted through analytical 


considerations, consist of integrating Eqs. 3a-3c and 4 of the paper (n = 1, N) 

Rx{n) = — A-1-811 sin^ Ry{n), 

(65a) 

= A-1-811 sin^ Rx{n) — ^Rz{n), 

(65b) 

Rz{n) = -1-8^sin^ ^Ry(n), 

(65c) 

and 

d„ -1- X + d-Rxin) + DRz{n)] sina„ = /3^(5a„, 

(66) 


where (5a„ = a„_i — 2 q;„ + a„+i and the parameters D, 0, Cp, /i. A, and x are defined in the paper as functions 
of the interaction matrix elements Vij {i,j = 1,2), the energy difference between the ground and the first excited 
level, El — Eq, and the discrete vector potential coupling parameter j3. The earlier system of ordinary differential 
equations is integrated in time with a standard fourth order Runge-Kutta algorithm with constant time-step At, 
typically 10“^. Such small time-steps, and even smaller sometimes, are required to conserve up to high accuracy 
the total and partial probabilities as the compound system of qubits and the electromagnetic vector potential evolve 
in time. In all the simulations, periodic boundary conditions are used. Due to the particular shape (Lorentzian) of 
the EM vector potential pulse and the population inversion pulse in the qubit subsystem, very large systems with 
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N = 2^^ = 8192 and N = 2^'^ = 16384 have been simulated to diminish as much as possible the effects from the ends 
(i.e., the interaction of the pulse tail with itself, as it could be in the case of periodic boundary conditions). In some 
cases, it is necessary to simulate even larger systems, with N = 50,000. In the figures presented in the paper and 
here only part of these arrays is shown, for clarity. In order to observe two-photon self-induced transparent (TPSIT) 
pulses an{t) and the induced quantum coherent population inversion pulses Rz{n\t), the following initial coditions 
are implemented: for the former, the analytically obtained solution resulting for the given set of parameters, while 
for the latter all the qubits are set to their ground, Eq, state. In terms of the Bloch vector variables Ri, i = x,y,z, 
that initial condition is specified as: 


R^,y{t = 0)=0, R,{t = 0) = -l, (67) 

for any n = I,..., iV. In that case, the above mentioned pulses exist for velocities less than the corresponding limiting 
velocity for TPSIT media, i.e.. 


9 k n k 

V <c = 13^- = 2P\. 

uj A 


( 68 ) 


The last equality is valid only when the two-photon resonance condition uj = A/2 is satisfied. In Eq. (20), k and u 


denote the wavenumber and frequency of the carrier wave of the electromagnetic vector potential; while the latter is 
determined by the two-photon resonance condition, the former may vary within an interval which is restricted by the 
qubit-parameter-dependent condition 


2x'(En+Eoo) < (Ei-Eo)2, (69) 

that guarantees the wavenumber k being real. 

The role of decoherence in TPSIT in superconducting quantum metamaterials like the one investigated here can be 
observed in Figures and in which the amount of decoherence, quantified by the factor 7 , takes the values 0, 
0.01 and 0.1, respectively. The decoherence factor 7 depends solely on the interaction matrix elements it is given 
by 


7 = 2 


Vn - Ido 


00 


VlQ 


(70) 


i.e., it is equal to zero for V/o = ^11 while it assumes non-zero values for Vn > Vqo- The Figures |4][^ show snapshots 
of Rz{n;t) and a„(t) at several instants from t = 0 to t = 168, which are separated by 14 time units, along with the 
corresponding analytical solutions. All these profiles but the first are shifted downwards to avoid overlapping, while 
only part of the array is shown for clarity. Note that time increases downwards. In Figure [4[ for 7 = 0, the a„(t) 
pulse is observed to excite a population inversion pulse i?z(n; t) in the qubit subsystem (Figure^). The amplitude of 
Rz{n]t) gradually increases until it attains its maximum value close to unity, while at the same time it propagates to 
the right with velocity v'. In the figure, that occurs for the first time at t ~ 70 time units; subsequently it evolves in 
time while it keeps its amplitude almost constant for at least 56 time units. After that, its amplitude starts decreasing 
until it is completely smeared (not shown). During the time interval in which the amplitude of the Rz(n;t) pulse is 
close to the predicted one (i.e., close to unity), the metamaterial is considered to be in the almost quantum coherent 
regime. Note that at about t = 84 a little bump starts to appear which grows to a little larger {’’probability bump”) 
which gets pinned at a particular site at n ~ 300. Note that another such bump appears at n = 0 due to the initial 
’’shock” of the qubit subsystem because of the sudden onset of the a„(f) pulse. A comparison of the numerical i?z(n; t) 
profiles with the analytical ones reveals that the velocity of propagation v’, which is the same as that of the numerical 
a„(t) pulse (Figure]^), is slightly larger than the predicted one v {v' > v). The corresponding profiles for very weak 
decoherence, in which the value of the decoherence factor 7 = 0.01, are shown in Figure 2. For that weak decoherence, 
there are only slight differences in comparison with Figure 1 which actually cannot be observed in the scale of the 
figures. 

In Figure]^ the decoherence factor has acquired a substantial value of 7 = 0.1, so that its effects are now clearly 
observed in the population inversion pulse Rz{n; t), while the vector potential pulse a„(t) does not seem to be affected. 
In Figure]^, while the i?z(n; t) is still excited by the vector potential pulse an{t), it has a very low amplitude compared 
with that of the analytically predicted form. Indeed, that amount of decoherence only slightly changes the analytical 
Rz{n;t) pulse and that change is hardly visible in the scale of the figure. Note that the speed of the i?z(n;t) is the 


even 


same as in the case without decoherence (which is also the same with that of the a„(t) pulse, in Figures 4]|6) 
the unwanted ’’probability bumps” appear at about the same locations with almost the same amplitude and shape 
independently of the amount of decoherence. 

Next, the possibility for the existence of propagating two-photon superradiant pulses in SCQMMs for large values 
of /3 is numerically explored (Figure 4). In order to obtain that figure, the qubit subsystem is initialized with all the 
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FIG. 4: Numerical validation of the analytical expressions for two-photon self-induced transparent (TPSIT) 
propagating pulses in superconducting quantum metamaterials (SCQMMs) without decoherence, a, Snapshots of 
the population inversion pulse Rz{n; t), excited by the induced quantum coherence in the qubit subsystem by the electromagnetic 
vector potential pulse, in the absence of decoherence (7 = 0 ); the pulse propagates to the right (time increases downwards), 
b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse an(t), that exhibits significant 
broadening while its amplitude decreases as time passes by; the numerical and analytical pulses are shown in blue and red 
color, respectively. Parameter values: x = 1/4.9, j3 = 6, Vqo = Vu = 1 (■y — 0), Voi = Vio = 0.7, Ei — Eq = 3, and v/c = 0.7. 


qubits in their excited state, while the vector potential pulse assumes its analytically predicted form for the selected 
parameter set. A very large array with N = 50,000 is simulated, and the initial position of the (center of) a„(t) pulse 
is at n = —18, 750. The subsequent evolution produces the snapshots presented in Figure 4, from t = 0 up to t = 300 
time units (time increases downwards). These snapshots are separated in time by 20 time units and they are displaced 
vertically to avoid overlapping. The population inversion Rz (n; t) and the electromagnetic vector potential pulse o„ (t) 
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FIG. 5: Numerical validation of the analytical expressions for two-photon self-induced transparent (TPSIT) 
propagating pulses in superconducting quantum metamaterials (SCQMMs) in the presence of weak deco¬ 
herence. a, Snapshots of the population inversion pulse excited by the induced quantum coherence in the qubit 

subsystem by the electromagnetic vector potential pulse, in the absence of decoherence (7 = 0 . 01 ); the pulse propagates to the 
right (time increases downwards), b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse 
a„{t), that exhibits signihcant broadening while its amplitude decreases as time passes by; the numerical and analytical pulses 
are shown in blue and red color, respectively. Parameter values: y = 1/4.9, /3 = 6 , Vqo = 0.998, Vn = 1 . 002 , Vbi = Vio = 0.7, 
El — Eo = 3, and v/c = 0.7. 


are shown in Figures and[^, respectively. In that simulation, the parameter /? is much larger than that used in 
Figure 3b of the paper, i.e., here /3 = 30. The other parameters, except the pulse speed, here v = 2.5 > c, are the 
same as those used in Figure 3b of the paper. The spatio-temporal evolution of both pulses reveals a rather complex 
pattern, that is discussed below. In Figure 4b, the a„(f) pulse breaks into several pulses of different amplitudes and 
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FIG. 6 : Numerical validation of the analytical expressions for two-photon self-induced transparent (TPSIT) 
propagating pulses in superconducting quantum metamaterials (SCQMMs) in the presence of substantial 
decoherence, a, Snapshots of the population inversion pulse Rz{n; t), excited by the induced quantum coherence in the qubit 
subsystem by the electromagnetic vector potential pulse, in the absence of decoherence (7 = 0 . 1 ); the pulse propagates to the 
right (time increases downwards), b, Snapshots for the corresponding evolution of the electromagnetic vector potential pulse 
a„{t), that exhibits signihcant broadening while its amplitude decreases as time passes by; the numerical and analytical pulses 
are shown in blue and red color, respectively. Parameter values: y = 1/4.9, /3 = 6 , Vqo = 0.98, Vn = 1 . 02 , Vqi = Vio = 0.7, 
El — Eo = 3, and v/c = 0.7. 


velocities; the most important are discernible along the green (c) and the purple (b) solid lines. In the same subfigure, 
the red (a) solid line indicates the analytically predicted pulse speed for two-photon superradiance in superconducting 
quantum metamaterials. The a„(t) pulses emerging from the initial one affect the qubit subsystem in a peculiar and 
complex way: First of all, a very well-defined Rz{n]t) pulse is observed which follows the analytical solution (in red) 
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FIG. 7: Numerical exploration of the propagation of two-photon superradiant (TPSRD) pulses in supercon¬ 
ducting quantum metamaterials (SCQMMs) for large values of /3. a, Snapshots of the population inversion Rz{n;t) 
induced in the qubit subsystem by the electromagnetic vector potential pulse, a„(t) (blue), in the absence of decoherence 
(7 = 0 ). The analytically predicted forms are shown in red color, b, Snapshots for the corresponding vector potential pulse 
a„{t) at the same time-instants. Parameter values: x = P = 30, Vbo = Vii = 1, Vbi = Flo = 0.8, Ei — Eq = 3, and 
w/c= 2.5 


for quite some time (about ~ 140 time units). After that time the pulse (whose last snapshot is numbered 1) slows 
down while it gets narrower while it is propagating together with the an{t) pulse along the purple (b) line. Note that 
around the initial position of the a„(t) pulse, the population inversion Rz{n;t) exhibits a region of strong variability 
within its extreme values ±1 which expands in the course of time. At t ~ 100 a second Rz{n;t) pulse is observed to 
emerge from the strong variability region, which is significantly narrower than the analytical solution. That second 
pulse ((whose last snapshot is numbered 2) propagates slower than the first one (1), along with the a„(<) pulse which 
can be observed along the green line (c) of Figure [^. 





































































